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Abstract: A three-dimensional unilateral contact problem for articular cartilage layers at¬ 
tached to subchondral bones shaped as elliptic paraboloids is considered in the framework 
of the biphasic cartilage model. The main novelty of the study is in accounting not only 
for the normal (vertical), but also for tangential vertical (horisontal) displacements of the 
contacting surfaces. Exact general relationships have been established between the contact 
approach and some integral characteristics of the contact pressure, including the contact 
force. Asymptotic representations for the contact pressure integral characteristics are ob¬ 
tained in terms of the contact approach and some integral characteristics of the contact zone. 
The main result is represented by the first-order approximation problem. 

1 Introduction 

Biomechanical contact problems involving transmission of forces across biological joints are 
of considerable practical interest (see, e.g. [21 El DU [13]) . Many analytical solutions to 
the problem of contact interaction of articular cartilage surfaces in joints are available. In 
particular, Ateshian et al. [8] obtained an asymptotic solution for the axisymmetric contact 
problem for two identical biphasic cartilage layers consisting of a solid phase and a fluid 
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phase and attached to two rigid impermeable spherical bones of equal radii. Later, Wu 
et al. [Hj extended this solution to a more general axisymmetric model by combining the 
assumption of the kinetic relationship from classical contact mechanics [123 wifi 1 the joint 
contact model [8] for the contact of two biphasic cartilage layers. An improved solution for 
the contact of two biphasic cartilage layers in the axisymmetric setting, which can be used 
for dynamic loading, was obtained by Wu et al. [ 15] . 

An asymptotic modeling approach to study the contact problem for biphasic cartilage 
layers has been performed by Argatov and Mishuris in a series of articles (see [I, ,5J 7]). In 
particular, it was shown [3] that accounting for the tangential displacements is important in 
the case of diseased cartilage where the measurement of indentation depth may differ even 
as much as 10% in comparison with the healthy case. In [5], the unilateral contact problem 
for articular cartilages bonded to subchondral bones with a contact zone in the shape of 
an arbitrary ellipse has been considered, and a closed form analytic solution was found. 
Exploiting this exact result, Argatov and Mishuris [7] have performed perturbation analysis 
of the contact problem with approximate geometry of the contact surfaces. Other analytic 
solutions for the contact problem were found using the viscoelastic cartilage model for elliptic 
contact zone in [6j. A new methodology for modeling articular tibio-femoral contact based 
on the developed asymptotic model of frictionless elliptical contact interaction between thin 
biphasic cartilage layers was presented in [2]. The mathematical model of articular contact 
was extended to the case of contact between arbitrary viscoelastic incompressible coating 
layers. 

In this study we extend results obtained in papers [3J and [5j by considering the influence 
of the tangential displacements on the contact problem for cartilage layers with the contact 
zone of elliptic shape based on the biphasic material model. Note that the perturbation 
method proposed in [7] could be one of the options for the analysis, however, the procedure 
is too complex to perform even a few asymptotic steps. Here, employing some technique and 
ideas from [3] and [5J, we propose another way to construct the asymptotics which utilises the 
assumption that the shape of the contact zone is an ellipse at the initial stage of deformation 
and can be regarded as a small perturbation of the ellipse at any other stage of deformation. 

The paper is organized as follows. The unilateral contact problem formulation and its 
linearization are presented in Section [21 where a special case of the contact configuration 
with one cartilage layer being plane and rigid is also considered in detail. In Section [3l 
we derive exact general relationships between the contact approach and some integral char¬ 
acteristics of the contact pressure, including the contact force. In Section 13.31 we obtain 
asymptotic representations for the contact pressure integral characteristics in terms of the 
contact approach and some integral characteristics of the contact zone. The zero-order and 
first-order asymptotic approximations for the solution to the contact problem are obtained in 
Sections 14.II and [4.21 respectively. Namely, the first-order approximation problem constitutes 
the main result of the present study. 
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2 Formulation of the contact problem 


We consider a frictionless contact between two thin linear biphasic cartilage layers firmly 
attached to rigid bones shaped like elliptic paraboloids. In the Cartesian co-ordinates 
(xi ,x 2 ,z) = (x, z) the equations for the two cartilage surfaces can be written in the form 
z = (-l) n $ (n) (x), n = 1,2, where 


$W( X ) 


x\ 


2 R\ 


(n) 


+ 


X% 


2 R^ ] 


( 2 . 1 ) 


with Ri'\ R 2 ^ being the curvature radii of the n-tli bone surface at its apex. 

In the undeformed state, the cartilage-bone systems occupy convex domains 0 ( x ) 

and 0 > <fA 2 ^(x), respectively. They are in the initial contact with the plane z = 0 at the 
origin of the co-ordinate system. 

We denote by Wi(x.,t), tu 2 (x, t) the local vertical displacements of the corresponding 
cartilage surfaces. Let also ui(x, t), u 2 (x, t) be the local horizontal (tangential) displacements 
of the corresponding surface of the cartilages. Finally, we denote by P(x, t) the contact 
pressure density. In this notation the equations for the cartilage surfaces can be written in 
the following form: 

z = Si(t) - (x + Ui(x, t)) + Wi(x, t), 


Z = — 62 (f) + <f> (2) (x + u 2 (x, t)) - iu 2 (x, t). 


( 2 . 2 ) 


Here, 5i, h 2 are some (positive) vertical displacements of the rigid bones. Note also that 
the vertical displacements w±, w 2 are positive, while the tangential displacements Ui, u 2 are 
directed outside of the contact zone. Denoting by <5*(t) = 5i(t) + 62(f) the contact approach 
of the bones, we get from (12.21) the following inequality: 


5*(£) + wi(x, t) + u> 2 (x, t) < $ (1) (x + ui(x, t)) + <h (2) (x + u 2 (x, t)) . (2.3) 


It was shown in [8j (see also 0) that vertical and the tangential displacements of each 
bone can be represented in the form 

w«(x, t) = { AP(x, t) + -^~ f AP(x, T)dr \ , n = 1, 2, (2.4) 

Un(x, t) — -VP(x, i), n = l,2. (2.5) 

Here e n = h n /a$ are dimensionless small parameters, hi, fi 2 mean the thicknesses of the 
cartilage layers, and a 0 denotes a characteristic measure of the contact zone (see the detailed 
description of the role of this parameter in [4]), H n = (A S) „ + 2 li s ,n)/Us,n are material pa¬ 
rameters of cartilages, where A Si „ and /i S;n represent the first Lame coefficient and the shear 
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modulus of the solid phase of the n-th cartilage tissue. Note that ui and U 2 in (12.51) do not 
necessarily coincide, they depend on both spatial variables X\, x 2 , and on the time variable 
t. 

Following [8], we introduce new spatial variables and time variable via formulas 




a 0 


3 = 1,2, 

ho 


where 


3// a ,lfci 3^i Sj 2^2 hs, 1 l^s, 2 

^ h\ h' 2 1 ^ A Sj i + 2[i Syl A Sj 2 + 2/i S]2 ’ 

a 0 is a characteristic measure of the contact zone, and k \, &2 are the cartilage’s permeabilities. 
In these variables we have the following relations on the contact area oj(t) encircled by the 
curve T(t) = doj{t)\ 


wi(x', F) + w 2 (x', f) = ( ) { AP(x', t') + x / AF(x', T r )dr' )> , (2.6) 


$< n) (x' + u n (x',t')) - $ (n) (x') - ^^V<f> (n) (x') ■ VP(x',f'), n = 1,2. 

2/is,n 


Further the equality in (12.31) , i.e., 

<5*(t') + Wi(x',F) + w 2 (x.',t') = $ (1) (x + Ui(x, t)) + <f> (2) (x + u 2 (x, t)) , 


(2.7) 


( 2 . 8 ) 


determines the contact area u>(t). 

Now we substitute (12.61) . (12.71) into (12.81) and obtain the governing equation relating the 
contact pressure with the vertical approach of the bones <5*(f) in the following form (from 
now on we keep the names of new unknown functions, e.g. <F(x) := <F(x'ao) etc.): 

t 

AP(x, t) + y J AP(x, r)dr = m ^<F(x) — <5*(t) — V<F(x) ■ VP(x, t) j . (2.9) 

o 


Here we have introduced the notation 


m = 


hj h 3 2 \ 

3 / r s ,i 3 113,2) 


( 2 . 10 ) 


$(x') = <F (1) (x) + $ (2) (x')- 

Thus, it follows from (12. ip and (12.lip that the functions $ and $ are given by 


( 2 . 11 ) 


$(x) = <h(xi,a: 2 ) = Ax\ + Bx\ 


( 2 . 12 ) 
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with 


1 


1 


1 


1 


A = 


2R[ l) 2 r[ 2) 


B = 


2 2 R {2) 


and 


$(x) = $>(x\,X 2 ) = Ax 2 + Bx 


(2.13) 


Note that the coefficients in A and B are positive dimensionless numbers, which are less 
than unit. 

Without loss of generality, one can assume that A > B. Then, Eq. (12.91) can be 
rewritten in an equivalent form, using all dimensionless parameters! 


AP £ (x,t) + y / AP e (x,r)dr = //(^(x) - 6 e (t) - eV$ 2 (x) • VP e (x,t)), (2.14) 


where the following notation has been introduced: 


'&j{*)=x[ + e 2 jX 2 , j = 

1,2, 6 e (t) = 

(2.15) 

H = Am , e\ = \/B/A , 

e 2 = \J B/A , £=j- 


It is important to note that 



X = 0( 1), 

He < x- 

(2.16) 


Discussion of the characteristic values of the introduced parameters is presented in Section 
[5] (see also (8, Tj). 

Since the solution of (12. 14(1 depends on the parameter e, it is customer to denote an 
unknown contact pressure by P = P £ in what follows. Note that the problem for e = 0 
coincides with that considered in [5], where an exact solution to this problem was found. 

Equation (12. 14)1 is the equation for determination of the contact pressure P e (x, i) > 0, 
x G c o e (t). In particular, in the case when the contact domain is represented by an ellipse 


UJ e {t) 


jx G M 2 : 


x 2 i P 2 (t,e)xl 1 
b\t,e) + b 2 (t, e) ~ /• 


(2.17) 


We supply Eq. (12.14)1 with the following boundary conditions: 

x G r(f), (2.18) 

x G T(f). (2.19) 

Note that in the axisymmetric case formula (12.141) coinsides with formula [31 (8)]. 


P e (x,f) = 0, 
dP 

w £ (-^) = o, 
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The equilibrium equation 

JJ P e (x.,t)dx. = F(t) 

U e {t) 


( 2 . 20 ) 


connects the external load F(t), unknown contact pressure P e (x, t), and unknown contact 
domain 0 J e (t). 


2.1 Special case of the contact configuration 

In order to check the content of formula (j2.9[) we consider here a special case, namely, we 
suppose that the lower part cartilage layer is plane and rigid (the same assumption was 
employed in [E]), it means that /i S;2 = oo and = oo, i.e., 

$ (1 ) =0, <h = <h (2) . 


In this case we have got the following equation for determination of the contact domain cn(f) 
in the form similar to (j2.9[) : 

t 

AP(x, t) + x f AP(x, r)dr = m ^$(x) — 5*(t) — V<h(x) ■ VP(x, t) j . (2-21) 


Here we will have 


3/^ 2 3/i s 

m = i x = 


hi ’ 


hi 


( 2 . 22 ) 


At the same time, small changes have to be made in the right-hand side of Eq. (I2.21[) as 
follows: 


r 2 ry.2 

*(*) = rrk + 


2 

( 2 ) 


2R ( ' 2) 2 R 2 

^ = hlapxj + h 2 2 apx\ 


2 ^R? 2 } ' 


Thus Eq. (12 .2 ljl can be rewritten as 


A?(x, t) + 


3/i s ,2^> 

hi 


7 AP{x ’ T)dT= 3 -f{^ 

n \ 1 


+ 


Xn 


2 R. 


( 2 ) 


s*(t) (2.23) 


3ao 


Xl d xi p{x,t) + ^d xl P{x,t) 

R-2 


R\ 


( 2 ) 


(2.24) 


It can be easily checked that in the axisymmetric case Eq. (12.23ft reduces to the governing 
differential equation obtained in [4]. 
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3 General relationships between the solution compo¬ 
nents 


3.1 Determination of the contact approach 

In our model we assume that the external load is non-decreasing. Thus, the contact domain 
is monotonically expanded, i.e. 


UJ £ (ti) C UJ £ (t 2 ), Vti < t 2 . (3.1) 

It is convenient to suppose also that the contact pressure is defined on the whole plane. For 
this we simply extend the density P e (x, t) by assuming that 

P e (x,t) = 0, Vx ^ u) e (t). (3.2) 

Integrating (j2. 14[) over contact domain co(t), we get 

t 


AF(x, t)dx + x JJ J AP e (x, r)drdx = 

Uj(t) Uj(t) 0 


= 1-1 JJ (^i(x) - 6 £ (t)) dx - £fi JJ VT 2 (x) ■ VP e (x,t)dx. 

Lj(t) Uj(t) 


(3.3) 


For simplicity of notation, we omit here (and everywhere in the next two sections) the 
subindex e in uj £ . 

From the monotonicity of the contact domain (13.ip and assumption (13.2(1 . it follows that 
the second integral on the left-hand side can be written in the form 

t t 



AF(x, r)dr(ix = J Jj AP e (x, r)dxdr. 

uj(t) 0 0 a j(t) 

Using the second Green’s formula 


(3.4) 


(w(x)Au(x) — v(x)Au(x)) dx = j 
u(t) r {t) 




with u = 1 and v = P £ (x,t ) we get the following relation in view of the boundary condition 

02T9]): 

r QP 

AP e (x, r)dx= / -^A(x, s)ds — 0, Vr < t. (3.6) 

w(t) r (t) 


Therefore, the both integrals on the left-hand side of (13.31) vanish 
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Further, we use the first Green’s formula 



(< pAil> + Vy ■ V^) dx = 


dip 

^ ds 


m 


(3.7) 


with ^(x) = T 2 (x) and <p(x) = P £ (x,t). In this case the integral on the right-hand side 
vanishes in view of (I2.18p . and we obtain the relation 


V'F 2 (x) • VP £ (x, t)dx = — JJ P £ (x, t)AT 2 (x)dx = — 2(1 + el)F(t), (3.8) 

Uj(t) ij(t) 


where we used the equilibrium equation (12.20)1 and the identity 

A'I' 2 (x) = 2(1 +e 2 ) 


(3.9) 


with e 2 being dehned in (12.151) . 

In what follows, it is convenient to have the following notation for the integrals of the 
product of k-th power of the function Ti and Z-th power of the function T 2 : 


4w(w) = // lf(x)*i,(x)dx >0, k,l = 0,1,2,. 


(3.10) 


In particular, ^4 0 ,o( w ) is fl ie area °f the contact domain. It is to remember that the constants 
Akj(u>) depend hnally on t, but we omitted this fact in the notation in order to avoid 
cumbersome expressions. Computations of A k) i(oj) for the elliptic domain (12.171) we included 
into Appendix (see Section 16.11) . 

Taking into account Eqs. (13.61) and (13.81) . we get 


m 


Aifi(uJ £ (t)) 2(1 + epe 

Ao fi (u e (t)) + A 0 ,oMZ)) 1 J 


(3.11) 


This formula allows us to compute the contact approach 5 e (t) as a function of the total 
external force F(t) and the main axes of the ellipse describing the shape of the contact zone, 
which in fact depends on time too. 


3.2 Some integral identity for the contact pressure 

In order to write out a more informative equation for the contact load, we use the following 
trick. We multiply both sides of (12.14)1 by the function v(x) = v k 2 (x) and integrate the 


















obtained equation over the contact domain uj(t) 


T 2 (x)AP e (x,f)dx + x 



u(t) 


T 2 (x)AP £ (x, r)drdx = 

o;(£) 0 

= A 4 jj 'h 2 (x)'I' 1 (x)dx -/i5 £ (t) Jj v h 2 (x)dx 

o;(t) cj(t) 

- ^ Jj ^2 (x)V$ 2 (x) • VP £ (xd)rfx. 

w(t) 


(3.12) 


Let us calculate the integrals in this relation by using Green’s formulas. For the first integral 
on the left-hand side we use formula (13.51) with u = T 2 , v = P £ and the boundary conditions 
(12.181) . (12.1 9p . Hence, we obtain 


'L 2 (x)A P e (x,f)dx = jj AT 2 (x)P e (x, t)dx. 

bj(t) Lj(t) 


Now taking into account (13.91) . we get 


jj 'f 2 (x)AP E (x, t)dx = 2(1 + e\)F(t). 

u(t) 


(3.13) 


For the second integral on the left-hand side, we apply the same approach, but interchange 
first the integrals over c o £ (t) and over r E (0, t) exploiting the load monotonicity. Therefore, 
we arrive at the equation 

t t t 



T 2 (x)AP £ (x, r)drdx 


T 2 (x)AP e (x, r)drdx = 2(1 + el) 


F{r)dr. (3.14) 


uj(t) 0 0 uj(t) 0 

For the first and second integrals on the right-hand side, we simply use the notation (13.10H , 
which gives 


Ti(x)T 2 (x)dx = Ai i(b) (3), 


T 2 (x)dx = H 0 ,i(6;/3). 


(3.15) 


Uj(t) Lj(t) 

Finally, for the third integral on the right-hand side, we make use of the following simple 
formula which follows immediately from the definition of Ty 


T 2 VT 2 
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Then we can apply Green’s formula (j3.7[) and the boundary conditions (12.181) . (12.19ft to find 
JJ 'h 2 (x)VT 2 (x) • VP £ (x,t)dx = -i jj AT 2 (x)P e (x,t)«2x. 


u(t) 


U}(t) 


By applying the second Green’s formula (13.5[) with u = P £ , v = T 2 , and the boundary 
conditions (I2.18p . (I2.19p . we represent this integral in the form 


*!(x)Vfi(x) • VP t (x,i)it = -i jj $5(x)AF ( .(x,«)dx. 


(3.16) 


Lj(t) Uj(t) 

This integral still contains the unknown density of contact pressure P £ (x, t). Let us define 

M^P £ {t) = JJ T^(x)AP £ (x,f)dx. (3.17) 

Oj(t) 

Now we rewrite the relation (13.121) by using the results for all integrals (I3.13p ~ 03.16p in 
the following form: 

2(1 + e 2 )/CP(t) = fiAi } i(u £ (t)) — ii8 £ (t)Ao ) i(uj £ (t)) + ^-At^P £ (t). (3.18) 


Here, we have introduced the Volterra operator K, as follows: 


t 

KF(t) = F{t) + x J F(r)dr. 

o 


(3.19) 


Note that the integral in the right-hand side of the equation (I3.18P allows to continue 
the same procedure to deliver an asymptotic estimate for this equation. 

We continue to proceed with Eq. (I3.18P on the next steps. 


3.3 Asymptotic estimates of the integral characteristics A i^P £ (t) 

Now we proceed to calculate the last integral in (13.181) . For this we multiply the governing 
integral equation (I2.14p by T 2 (x) U — 2) and integrate over contact domain cu(t): 


U}(t) 



^ 2 ( x )AP £ ( x,f)dx + X / ^ 2 ( x )AP £ ( x,r)drdx = 


j(t) o 


= h JJ ^ 2 (x)^i(x)dx-^ £ (f) JJ T J 2 (x)dx 

U}(t) uj(t) 

-tie [I T J 2 (x)VT 2 (x) • VP £ (x,t)dx. 


Uj(t) 
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(3.20) 




























(3.21) 


By using the same argument as on the previous step, we get 

P e (t) = nAij - fj,6 £ (t)A 0jj (a; (3) — /ie jj ^(x)V'h 2 (x) • VP e (x,t)dx. 

u(t) 

For the last integral we use the relations 

1 


*j(x)V*,(x) = —V« +1 (x) 
J + 1 


and 


jj V^ 2 +1 ( x ) • VP e (x, t)dx. = - jj A'I / 2 +1 (x)P e (x, t)dx. 

Uj(t) co(t) 

Therefore, the integral 

M®P e (t) = /PC" 1 U ld (u E (t)) - 5 e (t)Ao d (u e (t)) + -^KM^Peit) 


(3.22) 


has been obtained as a solution of the integral equation (13.211) . Here the inverse operator 
/CP 1 is defined by the formula 


K r 1 r(i) = Y(t) - x / Y(r)e- x{t - T) dr. 


(3.23) 


Performing the same computation, we obtain the following representation for the integral 
in the right-hand side of (13.181) : 


N 


2P 


-l 


{A 1J+1 (u,(t)) - 4(i)4u+iM0)} 

2e N r /r”xi* +! y(t). 


+ 


(N + 2)! 

Substituting this representation into Eq. (j3.18]L we finally get 

N j 

2(1 + e 2 2 )KF(t) = V Mu+i (".(<)) - 4 MA« + iK(t))} 


(3.24) 


or equivalently 


j=o 


+ 


N 


-N+l 


(N + 2)! 


fJ N+l ]C -N M (N+2)p e ^ 


(3.25) 


2(1 + e’)AC w+1 F(i) = g {A lJ+1 (u,,(t)) - 4WA>j + iM«))} 
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+ (A r + 2)! , ‘" +1 ' M( ' V+2> ' P * (t> ' (3 ’ 26) 

The latter relation allows us to determine the problem parameters asymptotically with any 
prescribed accuracy. 

Note that apart from the fact that the shapes of the contacting bones are elliptical 
paraboloids, no additional assumptions on the shape of the contact zone have been made. 
On the other hand, no proof was offered to show that the contact zone is approximately 
represented by an ellipse. This will be done later. 

Remark 1. For every t for which the contact pressure P e {t) is bounded and the contact 
region oj(t) belongs to a bounded domain, the remainder pv +2 )i H N+1 M. ( jV + 2 ) P £ (£) in formula 
(13.261) tends to zero as N —» oo. Thus, the series corresponding to the sum on the right 
hand-side of (13.261) is converging. 


4 Asymptotic solution to the contact problem 

4.1 Zero-order approximation 

First, we get solution of the problem for e = 0. In this case Eq. (I2.14p has the form 

t 

A?(°)(x, t) + x j AP^(x, r)dr = fi (\E f 1 (x) — , (4.1) 

o 

where l'i(x) is defined in (12.151) . Since we know from [5] that the contact zone is an ellipse 
at this stage of approximation we will have 

6. = < 4 ' 2 > 

-4o,o(k>o(£)) 


Using formula (14.21) and calculations presented in Section l64~l (see formula (16.61) ). one 
can find that 

7r jr 

A,o(^o(^)) = ~r-, di i0 (o;o(t)) = —^| (^o + e i) > (4-3) 

P 0 4Po 


and therefore 


S (0 \t) = 


bj (ft + £i) 

4 PI 


(4-4) 


Note that formulas (14.3j) and '11.11) contain two known constants e\ and e 2 defined in (12.151) 
and two still unknown functions bo(t ) and flo(t), which are the main semi-axis and the 
eccentricity of the ellipse 

„2 R2(+\ rr ,2 'i 

(4.5) 


H,(t) = ■■ X 6 » 2 : 4U + AU2< 1 


bo(t) b 2 0 (t ) 
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The leading terms in (13. 26ft imply (for N — 0) the following equation: 



2(1 + e 2 )JCF{t) — /iA 11 (o;o(t)) — /ih < '°'(t)A 0) i(a;o(t)). 

(4.6) 

Here. /C is the Volterra internal oDerator defined in (13.191). 

Analogously, using some results from Section 16.11 (see, in particular, formula flb.bDh we 
obtain 

-jrh 4 

Ao,i(^o(t)) — 3 (/3 0 + e 2 ) (4.7) 

4 Po 

and 

7 T h e 

Ai,i(o;o(t)) — 2 ^ 5 {3/5 0 + (ef + e 2 )(3 0 + , 

(4.8) 

and thus 

nb 6 

2(1 + e 2 )/CP(f) — {3/?o ( e i + e 2 )A) + . 

(4.9) 

To find the functions bo(t) and f3o(t) together with the pressure distribution 
contact zone, P^(x,f), we follow |5] and introduce a new unknown function 

over the 


t. 

p (0) (x,f) = P (0) (x,t) + y J P (0) (x,r)dr = /CP (0) (x,t). 

(4.10) 


o 


In the case of monotone external load, this function should satisfy the Poisson equation 
(following from (12 .Oft ) 


Ap (0) (x,t) = fi (^i(x) - h (0) (f|) , xGwo(f), (4.11) 

with the boundary conditions (12.181) . (12.191) . 

It is customary to rewrite this relation in the form 

Go(x,t) = 0, (4.12) 

where 

GoCM) = G 0 (b 0 , (3 0 , 5 q) (4-13) 

= Ap (0) (x,t) - /i (Ti(x) - h (0) (t)) , x e u 0 (t). (4.14) 

Bearing in mind that the function 4/i(x) is a quadratic polynomial (compare with (12.15j) ). 
it is natural to look for the solution of such problem in the form of a polynomial in x±,X 2 of 
the fourth degree, that is 

P W ( b o, A),Po,xp) = rj 0 {t) ^1 - p- - ) Qo(xi,x 2 ). (4.15) 
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Note that the term in the brackets vanishes on the boundary uq, and thus the condition 
(I2.18P is satisfied automatically. 

In Section I6T2T it has been shown that Qo is a polynomial of the second order having the 
form 

%(*,.*,)=( 1-f-^y). (4-16) 

so that 

p {0) (x 1 ,x 2 ]t) = 77o(t) ^1 - ^ ■ (4.17) 

Taken into account this representation we arrive at the following relations (see Section l6T3]k 


Vo(t) = 




4(1 + A 


b 2 

2 v 0) 


Vo (t) = 






2(6 + 2 f3 2 ) 4(3 + fil) ’ 


Vo (t) = 


i j K e i 


lK e \ 


2 (2/3q + 6/3q ) 4/5q (1 + 3/3q 


(4.18) 


(4.19) 


(4.20) 


This system allows us to determine the unknown functions b 0 (t ) and j3o(t). Indeed, 
eliminating from the last two equations, we get a bi-quadratic equation defining the value 
of the parameter An i.e., 

3/3 0 4 + (1 - eM - 3e? = 0. (4.21) 

By definition, /3 0 is a positive parameter, thus the unique positive solution of (I4.2ip has the 
form 


A) — 


[e\ - 1) + \Je\ + 34e^ + 1 


6 


(4.22) 


Note that at the zero-approximation the parameter A) does not depend on time. The other 
parameter, r]o(t), can be computed directly from (14.19j) or (14.2011 . if one knows the remaining 
constant bo(t). Moreover, taking into account (I4.18P and (11.11) . one can use an equivalent 
formula 

+Wo 2 + e i) 


Vo (t) = 


(4.23) 


1 6/?o(1 + AD ' 

In the same way, one can offer, in addition to (14.4p . two equivalent representations for 
the indentation parameter 


<$<°>(t) = 


i + A 2 




3 + $'“'-' /3o 2 (l + 3/3„ 2 ) 


(4.24) 


Finally, the major semi-axis bo of the ellipse loq is determined as follows: 

1/6 

96/3 0 5 (l + e|) 


bo(t) = 


F{t ) + x / F(r)dT 


V n (3f3o - /?o( e ? + e !) + 3efe|) 


(4.25) 
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Note that the parameters bo, Vo as well as the indentation, ho, depend on time t in contrast 
to the ellipse eccentricity /3 0 . 

Now, it remains only to find the pressure over the contact area. Using (14.10P and (I4.17[) . 
we get 

P W (bo,f3 0 ,rio,x 1 ,X2,t ) = /C _1 (vo{t)Qo( x i,x 2 ) 2 ) ■ (4.26) 

If (xi,x 2 ) belongs to the initial contact zone, i.e. 1 — -Prh- — > 0, then 


P {0 \x 1 , X 2 l t) = T) 0 (t) ^1 - - 


PiAY 
h l {t)J 


~ X / Vo (r) 1 - 


t .2 

X 1 


PlA x 2 


b 2 o(r) 6g(r) 


x 2 6 2 x 2 

If (xi,x 2 ) lies outside of the initial contact zone, i.e. 1 — y — < 0, then 


e~ x{t - T) dT. 

(4.27) 


U (xi,X2) 

The critical moment of time f* is determined by the formula 

b l(Q = A + Po x l- 

Using (14.251) . we get 


x; 


bl(T) 


a2 r 2 \ 2 

PoX 2 ' e-^-^dr. 


(4.28) 


F{t*) + X F(r)dr = 


I ' % - Me ' + ^ + 3e ‘ e b (*? + ^l) 3 . (4.29) 


96/3q 


1 + e. 


If the load is stepwise, we have F(t) = F () . Hence, we hnd that 

34 alt XI , „2\ , o„2„2 




/i7T 


VSPoXFo 

Note that in this case 


(3/3 0 A)( e l + e 2) + 3 e i e 2 ) , 2 o2 2\3 
n , 2 "T - Po X 2 ) 

i + e 2 


1 

x' 


fc) = 


96/3q (1 + e|)(l + yU 


2 „2\ *b. 


/z7t(3/5q - /?o( e i + e i) + 3e?e|) 


(4.30) 


(4.31) 


This hnishes the zero iteration step. Note that the results of this Section after changing 
the notation coincide with those obtained in [4j. 
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4.2 First-order approximation problem 

For the next steps we consider an appropriately deformed contact domain uj^\ defined as a 
perturbation of the zero-order one loq. Namely, we assume that it can be written in the form 

<4 1} = u;W(t) = |(xi,x 2 ) : Q 0 (x,t) + eQ 1 (x,t) > oj, (4.32) 

where unknown polynomials are taken in the forms 


Qo(x,t) = Q 0 (x,j3 1 ,b 1 ), 

Qi(x,t) = a 40 (t)xf + a 22 (t)x\x\ + a 04 {t)x 


(4.33) 

(4.34) 


Note that for e = 0 the solution form coincides with (14.51) . if one take b 4 = b 0j /3i = f3 0 . 
The idea behind such choice of the asymptotic anzatz is to satisfy the boundary condi¬ 
tions (I2.18P and (12.191) automatically. This will be archived by putting 


P. ll} = C 1 (rj il \t)(Qi 1 (x 1 ,X2,Pi(t),b 1 (t))+EQi(x,t))' 


(4.35) 


Now, when the boundary conditions are valid, we will satisfy the governing equation 
(12.9j) . Note that 

P £ (1) = Po + eP^Oie 2 ), (4.36) 


where pj = K,(Pj), j = 0,1, and 


= (D {t) I i _ 

1 bj(t) b\(t) 


Pi = 2 p {1 \t) ( 1 


x. 


PHf)x\ 


m m 

Substituting this representation into Eq. (12.91) . we obtain 


Qi( x i t). 


(4.37) 

(4.38) 


/C (A(?(°) + ePi + 0(e 2 ))) = p (tf i - 5™ - • (VP (0) + eVP (1) + 0{e 2 ))) , (4.39) 

where the parameter hg 1 ' ) is represented in the same form as Pj 1 ^, i.e., 

= So + + 0(e 2 ) = h (1) + 0(e 2 ). (4.40) 

We can write Eq. (14.391) with the accuracy to the terms of 0(e 2 ) as follows: 

A p (0) + eA Pl = p (Ti - <5 (1) - eVT 2 ■ VP (0) ) . (4.41) 


An extended variant of this equation can be written by using the definition of all com¬ 
ponents of the equation and by comparing coefficients at different powers of x i,x 2 , so that 
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Ar]^ 




4^( 1 ) 


3 + PI 

bt 


+ £( 6 d 4 o + CI 22 ) 


Arj^ 


r A 2 ( l + 3ft 2 ) 


&! 


+ £((3 22 + 6004 ) 


— p(l — 8£0 2 ,o), 

= /i(e 2 - 8 eel 6 2 , 2 ), 


24t 7 ( 1 > 


(040/^1 +022(1+^4)+ <304) —85/3(1+63)04,2, 


4 7 i( 1 ) 

— £ (<34o(15 + ^4) + 022) = 8£/30 4jO , 

4^(!) 

— £—^2 —(004(15/^4 + 1) + <322/^4) = 8 £/ 3 e 204 , 4 , 


where 


(4.42) 

(4.43) 

(4.44) 

(4.45) 

(4.46) 

(4.47) 

(4.48) 


02 fc, 2 i(<) = (r/ ( 1 ) 6 r 2fe ^) , k,l — 0 , 1 , 2 . 

In the system (I4.42IMI4.47I) we have 6 equations and 7 unknowns: <5 £ , bi(t), /3i(t), 

and 040 , 022,004 (coefficients of the polynomial Q 4 ). Therefore, we have to add an extra 
equation to the above system, namely 


6 U{t) = f + 

A),o(<^e(£)) A 0 ,o(o; £ (t)) 

where iq (f) can be represented in the form 


(4.49) 


uw = / / u (1) (x,()d 


x. 


„(i) 


We also make use of Eq. (13.261) written for this approximation step with the accuracy of 
0 (£ 2 ) in the form 

2(1 + e 2 2 )K 2 F{t) = g Mj.WkW {4 u+ 1 ( We (t)) - ^>(i)4„,, + 1 K(i))} • (4.50) 


Remark 2. Note that putting £ = 0, the system (I4.42MI4.47I) . (14.491) transforms to the 
previous case evaluated in the previous section. 

Remark 3. In the case when £ > 0, the system fl 1. 12D — (1 1. 17D . (14.491) has to be solved 
numerically. Note that the parameter £ in the last three equations (14.451) - (14.471) can be 
canceled. We left these multipliers here to explain the limiting case {e = 0). 
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5 Discussion and conclusion 


First of all, observe that at t — 0, the contact problem for biphasic layers reduces to that for 
elastic incompressible layers. The contact problem in the latter case were studied in a number 
of papers m El EDI Eg, however, without taking into account the tangential displacements. 

To solve the resulting problem 114.420 114.470 and (14.49H , we suggest the following iterative 
algorithm: 

• Taking £ = 0, we have computed all values 77 , b,/3,6 = rj 0 , b 0 , /3 0 , So from the zero-order 
approximation. 

• Having them we can compute the quantity 02 k,2i(t) from (14. 481) . 

• Then, from the system of three equations (I4.45IMI4.47I) we compute the constants 
CI 40 , 022 , o 04 assuming the values of rj, b, f3 as above. 

• Finally from the system of four equations (I4.42p - fl4.44p and il 1.1 91) considering the right- 
hand side known (computed by the values know from the previous computations), we 
found new values rj,b,/3,5 and compare them with the previous computations. If the 
required accuracy has achieved we stop the computation, if not we are going to the 
second step of this iterative procedure. 

We note that formulas (12.4p and (12. 5p for the vertical and tangential displacements 
contain different powers of parameters e, namely, e 2 and e, respectively. Note also that our 
analysis (with the values of another parameters taken into account) shows, that the role 
of these magnitudes (vertical and tangential displacements) is quite opposite. In the final 
equation (see (12. 14ft ) the leading terms, corresponding to the vertical displacement, contain 
the zero power of the new small parameter £, but the leading terms, corresponding to the 
tangential displacements, contain the first power of e. 
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6 Appendix 


6.1 Calculation of the constants Aj~i 

Here, we compute the values of the constants 


Ak,i(b ; P) 


^i(x)^2(x)rfx > 0, k,l — 0,1,2,... 


U}(t) 


First of all, we note that an unknown contact domain ou(t) is of the same type as sections of 
the initial gap elliptical paraboloid, i.e., it is an ellipse coaxial to the ellipse 


co(t) = u £ (t) 


jx G M 2 : 


x i , xfflfae) < A 

b 2 (t-s) + b\t-e) ~ /■ 


In order to avoid long formulas, we use the short notation for oj(t), writing all parameters 
without variables they depend on, i.e., 


u(t) 


. x 2 (5 2 


x e 


r 2 

fi + 

b 2 b 2 


< 



Performing the standard change of variables 

b 

x\ = br cos 9, X2 = — sin 9 , 


we represent the integral for Akj(t) i n the form 

1 2ir 


Ak,i(p'iP) — 



b 2 r 2 cos 2 9 + ^-^-r 2 sin 2 9 I ( b 2 r 2 cos 2 9 + “—^Ar 2 sin 2 9 I '—rdrdO 


b 2 ej 


,2 ■ 2 


0 0 


jj2k+2l+2 


P 2 

2t r 


2„2 „„„2 , 


b 2 e 2 2 2 . 2 X b 2 


P 2 


r 


2fc+2i+l 


dr I^a 

n i =0 


k] e l -„2 i 


^i\(k-i)\p 


2 i 


sin 2 * 9 cos 2k ~ 2i 9 


x 


V . uj n TTT ^ sin 2j 0cos 2/ - 2j ddd 


fr2k+2l+2 


(2 k + 21 + 2)/3 


2 77 ^ 


X 


V sin 2j 0cos 2/ - 2j ddd. 

P j\(l - ])\P 2 > 


P 


( 6 . 1 ) 


( 6 . 2 ) 
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Since the trigonometric functions are presented here only in even powers, then the last 
integration can be performed over the interval [0, 7 t/ 2] as follows: 


Ak,i(b;/3) — 


4f> 2 ‘+ 2 ‘+ 2 ?* k\ ef . 2ia 2k _ 2ia 

> —-- sill 6 cos 6 

Z-/ 7*1 


( 2 k + 2/ + 2)f3 J ^ i\(k-i)\p 2i 

o ?_0 


i 


x 


Eu 


n ej 

^p,k+2l+2 


sin 


2j Q COS 21 - 2 ! Q d Q 


2 i 


E k\ e 1 

, =o /!(/.- - i)! W 2 ‘ 


x 




/! 


7r/2 


j=z M - i) } -P 2 * J 

J- u o 


sin 2i+2j 0 cos 2fc - 2 * +2/ - 2j QdQ. 


(6.3) 


The integrals in (16.3j) are calculated by using formulas 

tt/2 

J sin 2p 0 cos 29 QdQ = ^ 
o 


2p fl ™2 8 ^_ir( P +i/2)r( 9 + i/2) 


r(p + g + i) 


, p, 7 > 0, 


(6.4) 


and Legendre’s duplication formula for the Gamma-function 


_. / x V2nT(2n) 

r(n + 1/2) = n G N, 


(6.5) 


2 2n- 1 /2r( n )' 

as well as the relation T(n + 1) = n\. Finally, we arrive at the following representation of 
A kt i = A k i(b ; /3) valid for all fc,l6N 0 = NU {0}: 

s2k+2l k 


A k ,i = 


27Tb 2 ( 6 / 2 )" 


/3(2k + 2l + 2)(k + /)! ^ i!(Jfe - *)! /3 2i 


Ett 


k\ e?* 


( 6 . 6 ) 


x 


E. 


l\ e 2 j ( 2 i + 2 j)\( 2 k- 2 i + 2 l ~ 2 j)\ 


“ j!(^ - i) ! /^ 2i (* + j)K k ~ i + 1 ~ i) ! 


6.2 Computation of the polynomial Qq 

In order to determine the coefficients of the polynomial 

Qo(xi,x 2 ) = 1 + 7i,c+i + ?o, 1^2 + g 2 ,o^? + qi, 1 X 1 X 2 + qo, 2 X 2 2 , 

we need to compute the normal derivative of the unknown functions p ^ 04.151) along the 
elliptic boundary V : 

dpi 0 ) 


dn 


r = Vp(»> ■ + r = !,„(() (-^1 - 2 -A d) Q„| r = 0. (6.7) 
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Here we take into account the fact that, since the contact domain is an ellipse (14.50 . the 
tangential and normal vectors to the boundary T = d£l are given by 


= (~P o x 2,^i) , "rt = (x 1: ^x 2 ) ■ (6.8) 

Then, to satisfy the boundary condition (12. 19jl the following equation should be valid: 


Qo|r — 0. 


This, in turn, is equivalent to the representation 


Qo(x 1 ,x 2 ) 


( xj Plx%\ 

V bl )■ 


(6.9) 


( 6 . 10 ) 


6.3 Evaluation of the ellipse parameters 

Since 


p (0) (x,t) = p {0) (x u x 2 ,t) = 7 l0 (t) ^1-^j-- 


Plx 2 


b 2 

°o 


we have 


<9p(°) 

dx 

d 2 p(°) 


0 ) / r 2 


PqX 2 

bl 


dx\ 


= 2//0 


b 2 


rr 2 «2 2 

X 1 Po x 2 

h 2 7,2 

°0 u 0 


+ 


2 x i 

°o 


2xi 2xi 
7,2 

u 0 u 0 


Therefore, by straightforward computations, we find that 


<9 2 p(°) 

dx 2 


= 2r/ 0 


2 ,6a:? 2/3 2 x 2 


6$ 


dp (0) _ 9 A _ x\ _ filx\ 
dx 2 ' l0 V bl b 2 




2 ffix 2 

bl 


d 2 p^°) 
dxl 


Thus, we obtain 


= 2770 


d 2 p(°) 


2A 


68 


1 _ 

ft 2 

°o 


nxf 2 \ 2f3$x 2 2f3$x 2 


b 2 
°o 


+ 


= 2770 (t) 


2/3 2 , 2/3gxf , 6/^1 


ft 2 

°o 


ft 2 

°0 




+ 


T, 4 

°0 


+ 


b 4 o 


Substituting (j6.13j) and (I6.15P into the main equation 

Go(b 0 ,/3o,S 0 ) = Ap (0) (x,t) - p (Ti(x) - 5 {0 \t)) = 0, 


( 6 . 11 ) 

( 6 . 12 ) 


(6.13) 

(6.14) 


(6.15) 


(6.16) 


where 


Go — 2?7q (t) 


l+ft 2 , (6 + 2/3A 2 76/3*+ 2/3 0 2 \ 

~ + v - ^- ) 1 + \~^r~) 2 


/i ( v ki(a;i,x 2 ) - <5 (0) (t)) , 
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and taking into account that 


'ki(x) = ^ 1 (x l ,x 2 ) = xj + e\x\, 

one concludes that the expression for G 0 is represented by a second order polynomial with 
respect to the independent variables X\ and X 2 in the following form: 

Go(b 0 ,p 0 ,Vo,8 {0) ) = qo(b 0 ,Po,Vo,8 {0) ) + Po,Vo)xj + q 2 {b 0 , Po,Vo)xl- (6.17) 

Here the coefficients are defined as follows: 

q 0 (b 0 ,f3 0 ,7] 0 ,5 {0) ) = ^§(1 + /3q) - <5 (0) , (6.18) 

qi(bo,Po,Vo) = ^x( 3 + /^o) ~ V, (6H9) 

°o 

92(^0) P 0 ; Vo) = (1 + 3^o) — V e i- (6.20) 

°0 


6.4 Auxiliary computation 

Taking into account (14.37j) . we can represent p o ( x )^0 in the form 


Po( x ,t) = V (1) (t ) ( 1 - ^ 


4^,4 
2 


bj 


2/3jxjx 2 2 xj Pfx: 

bf bj hf 


( 6 . 21 ) 


Hence, applying the Laplace equation, we get 

Ap 0 (x,f) = V {1 \t) ^-^(l + /?i)+x^(3 + /3 1 2 )+^^(1 + 3^)^) . (6.22) 

Next, by using representation (14.381) . we can write p 4 (x, t) in the form 

/ ,\ o (iL-a f 4 1 2 2 1 4 a -ioxj <222X1X2 ao 4 xfxi . 

Pi (x,t) = >(t) I + a 22 x 1 x 2 + a m x 2 -^-^-^— ( 6 - 23 ) 

_ aioPfxjxl _ a 2 2 Pjxjxj _ Oo4^|\ 
bj bj bj J ' 

Therefore, we obtain 


Api(x, t) = 2 ?7 (1) (t) ^(12a 40 + 2a 22 )a:f + (2a 2 2 + 12a 04 )^2 

12/3{a 40 + 12022(1 + /3 2 ) + 12a 04 2 2 
^2 

o 4 o(30 + 2/3 4 ) + 2022 4 2a22/3{ + Oo 4 (2 + 30/3 2 ) 

bj Xl bj x y ■ 


(6.24) 


(6.25) 
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We also use the following representations: 


^j{x) = x 1 + e j x 2 , j = 1 , 2 . 
Thus, applying the gradient operator, we simply get 

VT 2 (x) = (2x 1 ,2e 2 2 x 2 ) 

and 

VP 0 (x,t) = (/C- 1 Vp 0 (x,-)) (*)■ 

It yields the following representation: 


VT 2 (x) ■ VPq(x, t) - -8 /c: 


—1 


(!) f 1 _ _ 




X, 


+ 


e^i^s 


b\ J \bj bj 


(t) 


= -8xj 1C 


-' — l 


t ]W 

w 


(t) - 8 e 2 x 2 ( 1C 1 ) (*) + 8 - x 't 1 ^^4 - ) ) (t 


+8{l + e 2 2 )x{x 2 2 ( 1C 1 (^~Jr 


(t) + 8e\x A 2 ( 1C 1 (^-Jr 


(t) 


=: -8xl6 2 ,o (t) - 8elx 2 2 6 2)2 {t) + 8x A 9^ 0 (t) + 8(1 + e 2 )xfx 2 9 4 , 2 (t) + 8 e 2 2 x 2 9 4A (t). 
Here we have introduced the notation 


<k k ,2, = (1C-' W%T 2k ff)) ((). 

Combining the above results we obtain the system of equations fj4.42|) - fj4.47jl . 


(6.26) 
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